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Abstract 

Kalman  filters  and  ARIMA  models  provide  optimum  control  and  evaluation  tech¬ 
niques  (in  a  minimum  squared  error  sense)  for  clocks  and  precision  oscillators.  Typ¬ 
ically,  before  the  models  can  be  used,  an  analysis  of  data  provides  estimates  of  the 
model  parameters  (e.g.,  the  phi’s  and  theta’s  for  an  ARIMA  model).  These  model 
parameters  are  often  evaluated  in  a  batch  mode  on  a  computer  after  a  large  amount 
of  data  is  obtained. 

An  alternative  approach  is  to  devise  an  adaptive  algorithm  which  “learns”  the 
important  parameters  while  the  device  is  being  used  and  up-dates  the  parameters 
recursively.  Clearly,  one  must  give  up  some  amount  of  precision  if  one  deviates  even 
slightly  from  the  truly  optimum  techniques,  but,  as  this  study  shows,  the  costs  in 
performance  are  not  large  at  all.  If  one  chooses  the  best  sampling  intervals,  the  loss 
in  precision  can  be  negligible. 

The  physical  models  used  in  this  paper  are  based  on  the  assumption  of  a  com¬ 
bination  of  white  PM,  white  FM,  random  walk  FM,  and  linear  frequency  drift.  In 
ARIMA  models,  this  is  equivalent  to  an  ARIMA(0,2,2)  with  a  non-zero  average  sec¬ 
ond  difference.  Using  simulation  techniques,  this  paper  compares  real-time  estimation 
techniques  with  the  conventional  batch  mode.  The  criterion  for  judging  performance 
is  to  compare  the  mean  square  errors  of  prediction  between  the  batch  mode  and  the 
recursive  mode  of  parameter  estimation  operating  on  the  same  data  sets. 


INTRODUCTION 

Before  working  directly  on  the  ARIMA(0,2,2)  models  [1],  it  is  of  value  to  establish  a  few  important 
relations.  An  ARIMA(l,0,0)  model  is  often  referred  to  as  an  exponential  filter  since  its  impulse 
response  function  is  an  exponential.  That  is: 


—  4>Xn~\  4-  an 


(1) 


=  p 


for  n  =  1,  2,  ... 


Where  the  input  to  the  filter,  an,  is  taken  to  be  ao  —  1  and  an  =  0  for  n  >  0  (the  unit  impulse). 
This  filter  is  the  digital  equivalent  of  a  simple  R  -  C  low-pass  filter  with  <j>  —  exp(-t/RC)  where 
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t  is  the  time  interval  between  steps  in  the  counting  index,  n,  and  R  and  C  are  the  resistance  and 
capacitance  values  in  the  analog  filter.  In  frequency  domain,  the  filter  transfer  function  is  unity  for 
low  frequencies  and  drops  off  at  -6dB/octave  from  the  cut-off  frequency,  /„  ,  given  by: 

h  =  1/(27 vRC)  (2) 

Given  the  cut-off  frequency,  /„,  and  the  sample  time  interval,  t,  the  ARIMA  coefficient,  <j>,  can  be 
calculated  with  the  equation: 


<f>  =  exp(-2ic  fet0)  =  exp(-tofr)  (3) 

Thus,  for  a  given  physical  system,  the  ARIMA  parameter  is  dependent  on  the  data  sampling  rate. 

The  inverse  function  for  the  ARIMA(1,0,0)  is  just  an  ARIMA(0,0,1)  which  can  be  seen  by  solving 
Eq.  1  for  a„.  Such  a  filter  would  be  constant  for  the  low  frequency  end  changing  to  an  increase  of 
gain  by  +6dB/octave  above  the  same  cut-off  frequency  derived  above. 

The  ARIMA  models  considered  here  are  models  for  the  phase  of  the  clock  comparisons  (not  their 
instantaneous  frequencies).  We  can  begin  with  a  physical  model  of  a  clock  with  pure  White  FM. 
Since  the  ARIMA  model  is  for  phase,  this  first  example  would  be  given  by  an  ARIMA(0,1,0)  model, 
a  random  walk  of  phase.  The  power  spectral  density  (of  phase)  for  this  noise  is  that  of  a  constant 
decrease  of  -6dB/octave  everywhere  below  the  Nyquist  frequency,  l/(2to). 

We  can  add  white  noise  modulation  to  this  model  by  going  to  an  ARIMA(0,1,1)  model  where  the 
moving  average  parameter,  0,  is  computed  from  Eqs.  2  and  3,  above,  for  r  =  RC.  In  actuality  this 
is  not  an  adequate  model  for  many  real  clocks.  One  often  encounters  more  low  frequency  divergent 
noises  than  the  ARIMA(0,1,1)  which  require  an  additional  integration:  that  is,  one  needs  either  an 
ARIMA(0,2,1)  or  an  ARIMA(0,2,2).  Physically,  the  ARIMA(0,2,l)  model  corresponds  to  a  super¬ 
position  of  white  FM  and  random  walk  FM.  Again  the  transition  between  the  two  noise  regimes  is 
accomplished  by  using  Eqs.  2  and  3,  above. 

If  one  now  adds  white  PM  one  must  go  to  the  ARIMA(0,2,2)  model  (see  Fig.  1),  which  has  two 
break  points  corresponding  to  the  transitions  between  random  walk  and  white  FM  and  between  white 
FM  and  white  PM.  The  equations  above  allow  one  to  calculate  the  two  parameters  separately,  say, 
9[  and  0'2,  for  theta-values  of  two  cascaded,  MA  filters.  These  two  filters  can  be  combined  in  one 
MA(2)  filter  whose  theta-parameters  must  be  combined  as  factors  to  realize  the  correct  MA  filter. 
This  combination  is  obtained  as  follows: 

(1  -  0\B)(l  -  0'2B)  =  (1  -  0\B  -  02B2)  (4) 

which  yields  0\  =  (9\  +  0'2)  and  $2  =  -(^1^2)  where  B  is  the  index  lowering  operator  [l]  and  0\  and 
$2  are  calculated  using  Eqs.  2  and  3.  A  linear  frequency  drift  can  also  be  important. 

ADAPTIVE  APPROACH  TO  TESTING  ARIMA  PARAMETERS 

There  are  many  methods  of  estimating  parameters  —  for  example,  just  a  guess  is  one  means.  Of 
importance  are  issues  such  as  bias,  confidence  intervals,  efficiency,  and  likelihood.  While  it  is  easy 
to  present  various  estimating  procedures  that  “work”  it  is  often  difficult  to  evaluate  how  well  they 
perform.  This  section  develops  theoretically  the  effects  of  errors  in  the  ARIMA  parameters. 

Using  simulation  techniques,  these  errors  are  evaluated  in  "real-time”  and  compared  to  the  conven¬ 
tional  batch  method,  after  the  fact.  The  theoretical  consequences  of  parameter  errors  are  surprisingly 
mild.  That  is,  many  models  are  robust  in  regard  to  fairly  poor  estimates  and  they  can  give  surprisingly 
good  results.  Still  there  are  regions  of  operation  where  problems  can  arise.  For  example,  taking  data 
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very  frequently  does  not  improve  one’s  knowledge  of  basically  long-term  performance.  Knowledge  of 
the  long-term  performance  simply  requires  long-term  data. 

The  adaptive  (real-time)  estimation  is  based  on  the  fact  that  an  ARIMA(0,0,1)  has  an  autocorre¬ 
lation  function  given  by: 


jR(n)  = 


I 


-5/(1  +  52) 

0 


for  n  —  0 
for  n  =  1 
otherwise 


(5) 


Figure  2  presents  a  block  diagram  of  an  adaptive  algorithm  which  effectively  servos  the  theta- 
estimate  (denoted  “<£” )  to  the  “true”  theta-value  and  also  provides  a  current  estimate  of  the  variance 
of  the  residuals  and  a  drift  estimate.  Basically  the  algorithm  computes  contributions  to  the  first 
autocorrelation  coefficient  which  in  turn  adjusts  the  estimated  theta-value,  effectively  driving  the  first 
autocorrelation  coefficient  of  the  residuals  to  zero.  The  first  autocorrelation  coefficient  of  the  residuals 
(given  by  Eq.  8,  below)  is  proportional  to  the  difference  —  (see  Ref.  [l])  giving  both  direction  and 
value  for  the  servo.  Figure  3  depicts  a  simulation  of  the  servo  performance.  Although  the  servo 
“works”  we  need  to  compare  its  performance  with  a  conventional  batch  approach  for  estimation. 


THE  ARIMA(0,2,l)  MODEL 

Figure  4  is  a  diagram  of  the  forecast  system  for  an  ARIMA(0,2,l)  noise  model.  In  this  case,  theta  is 
the  “real”  parameter  which  is  unknown,  but  is  estimated  with  the  value  of  phi.  If  phi  were  to  equal 
theta,  then  the  system  would  provide  an  optimal  forecast  of  Xn.  Since  there  will  always  be  some  error 
in  the  estimate  of  theta,  it  is  of  value  to  explore  the  consequences  of  such  an  error.  The  following  is 
a  detailed  description  of  the  estimation  process  and  the  evaluation  of  the  errors  caused  by  an  error  in 
phi. 

Following  Fig.  4,  the  data,  Xn,  are  the  only  observables  from  the  clock  comparisons  and  the  model 
is  white  noise  FM  and  random  walk  FM:  That  is,  an  ARIMA(0,2,1)  model.  I  explicitly  assume  that 
the  model  is  a  good  model,  but  theta  is  unknown.  Theta  can  be  estimated  by  the  methods  given  in 
Ref.  [1],  or  by  recursive  filters  developed  here. 

The  output  data,  Xn,  are  filtered  with  an  “inverse”  filter  and  the  residuals,  Wn,  are  obtained.  In 
the  Box  and  Jenkins  method,  phi  would  be  adjusted  to  give  minimum  variance  to  Wn.  This,  however, 
is  accomplished  in  a  batch  mode  after  the  fact  and  not  in  real-time.  Regardless  of  how  theta  is 
estimated,  phi  is  used  in  the  forecaster  as  shown  in  Fig.  4:  Indeed,  there  is  no  other  value  to  use. 
With  Wn  as  input  to  the  forecaster  (the  switch  in  Fig.  4  in  the  “up”  position),  the  system  is  allowed 
to  run  for  a  time  to  let  all  transients  die  out.  At  this  point  the  estimated  output,  Xn  =  Xn  exactly 
since  the  inverse  filter  is  the  exact  inverse  of  the  forecast  filter. 

At  some  point  in  time,  t  =  n  +  1,  Wn  is  no  longer  available  as  input  to  the  forecast  filter  and 
its  input  is  set  to  zero  (switch  in  the  “down”  position).  For  the  first  forecast  value  Xn  +  1,  the  two 
previous  values,  Xn  and  Xn~i  are  available  to  use  in  the  forecast  as  shown.  The  error  of  the  forecast 
(after  the  fact),  can  be  found  by  subtraction  of  the  two  equations: 

■An+i  —  2An  _An_i  +  flfj-t-i  —  5a„ 

Xn+l  =  2  Xn-Xn.l~9Wn _ 

S\  —  —  An-t-l  —  ar»+l  — 


(6) 
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where  81  is  the  error  in  the  first  forecast. 

Following  Box  and  Jenkins  [l],  W„  can  be  expressed  as  an  infinite  series  which  incorporates  the 
psi-weights  and  the  uncorrelated  innovations,  an: 


Wn  =  E*<a»~i  (7) 

0 

The  psi-weights  can  be  expressed  in  terms  of  phi,  theta,  and  the  innovations,  an,  as  shown  below 
by  requiring  Wn  to  satisfy  the  equation  (see  Fig.  4): 


Wn  -  <t>Wn_ !  =  an  -  0a„- 1 


for  all  n,  the  result  is  given  by: 


-U 


for  i  =  0 

(<f>  -  for  i  =  l,  2,  3,  ••• 

We  can  now  evaluate  the  expected  square  of  the  first  forecast  error  as: 


(8) 

(9) 


e[«!]  =  <7![i  +  !  (io) 

For  phi  equal  to  theta  we  obtain  the  classical  result  that  the  variance  of  the  first  forecast  error  is 
just  the  variance  of  the  innovations. 

We  can  repeat  this  calculation  for  t  =  n  +  2  by  using  the  forecast  value  Xn+\  in  place  of  the 
(unknown)  value,  Xn+i.  Similarly,  for  t  =  n  +  3  and  so  forth.  The  result  can  be  summarized  in  the 
following  formula  for  the  mean  square  time  interval  error  for  M  lags  in  the  future: 

E\Sjg\  =  Ma\  {l  +  (1  -  6)(M  -  1)[1  +  (1  -  0)(2 M  -  l)/6]  +  M(<f>  -  0)2/(l  -  ^2)}  (11) 

This  formula  has  been  verified  using  simulation  techniques. 

There  are  two  points  to  make  in  regard  to  this  relation:  (l)  for  phi  equal  to  theta  the  result 
is  identical  to  the  classical  results  as  it  should  be,  and  (2)  sis  M  gets  larger,  the  variance  grows  as 
M-cubed  but  the  term  proportional  to  (<f>  —  0)  grows  as  M-squared.  That  is,  for  sufficiently  large  M, 
the  errors  of  parameter  estimation  become  unimportant.  Figure  5  shows  the  regions  of  forecast  errors 
(l)  primarily  due  to  conventional  analyses  and  (2)  those  due  primarily  to  an  error  in  the  estimate  of 
theta,  i.e.,  phi. 

Figure  5  shows  clearly  that  problems  develop  near  theta=l  and  near  phi=l  with  phi  and  theta  not 
near  each  other  in  value.  The  problems  near  theta  equal  to  1  can  be  reduced  by  having  longer  term 
data.  Having  more  frequent  data  doesn’t  help. 


COMPARISONS  BETWEEN  ADAPTIVE  AND 
CONVENTIONAL  PARAMETER  ESTIMATION 

Given  an  ARIMA  model  corresponding  to  a  physical  model  we  can  simulate  a  noise  sample  and  treat 
the  data  as  if  it  were  real-time  data  and  estimate  the  parameters  recursively.  The  same  data  can  be 
treated  in  batch  mode  and  find  those  estimates  which  minimize  the  sum  of  the  squares  of  the  residuals. 
Table  1  summarizes  the  results  of  the  estimation  process.  The  program  generated  100  noise  samples 
each  of  200  data  points  in  duration.  Each  sample  noise  was  processed  through  the  adaptive  estimation 
procedure  developed  here  and  the  conventional  Box  and  Jenkins  [l]  treatment.  Of  course,  the  “true” 
values  of  the  parameters  are  also  known  since  this  is  only  a  simulation. 
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TABLE  1  ARIMA(0,2,1) 
Theta  =  .9049  Sig-A  =1.1051 

Quantity 

Conventional 

Method 

Adaptive 

Method 

Phi  (Actuals. 9049) 

.8983 

.8577 

Std.  dev.  of  mean 

.0074 

Bias  rel.  Actual 

-.0066 

-.0471 

T-ratio 

-1.58 

-6.38 

Sig-A  (actual=1.105l) 

1.1161 

Std.  dev.  of  mean 

.0119 

Bias  rel.  Actual 

.0111 

-.0674 

T-ratio 

.93 

-4.63 

As  expected,  the  conventional  estimates  are  more  accurate  and  more  precise  than  the  adaptive 
methods  developed  here.  Still,  with  reference  to  Fig.  5,  even  fairly  large  errors  in  phi  relative  to 
theta  are  soon  covered  by  the  conventional  errors.  Table  1  shows  that  for  theta  =  .9049,  the  adaptive 
method  produces  statistically  significant  biases  (T-ratios  of -4.63  and  -6.38). 

ARIMA(0,2,2)  NOISE  ANALYSIS 

To  estimate  the  theta-value  the  servo  was  based  on  the  fact  that  an  ARIMA(0,0,1)  model  has  only 
the  zeroth  and  first  autocorrelation  coefficients  non-zero.  An  ARIMA(0,0,2)  model  has  an  additional 
non-zero  coefficient  at  lag  2,  which  is  strongly  dependent  on  the  second  MA  coefficient.  The  new 
servo  has  a  separate  loop  which  takes  samples  of  the  (lag  n)  (lag  n-2)  product  and  adjusts  <j> 2  to  null 
the  average  similarly  to  the  loop  shown  in  Fig.  2,  for  ^1. 

Similarly  to  Fig.  3,  Figure  6  depicts  the  transient  response  of  the  ARIMA(0,2,2)  adaptive  estima¬ 
tion.  As  noted  above,  the  important  performance  is  that  of  forecast  errors  not  the  intermediate  values 
of  $1  and  $2  shown  in  Fig.  6.  Still  it  does  show  how  the  parameters  stabilize. 

As  noted  above,  the  means  of  comparing  algorithms  is  to  compare  the  forecast  errors  for  similar 
data  situations.  The  Box  and  Jenkins  method  can  be  used  to  estimate  the  forecast  errors  using  the 
“psi-weights” ,  as  in  equations  8  and  9,  above.  The  model  for  consideration  now  is  as  ARIMA(0,2,2). 
The  theoretical  forecast  errors  for  the  model  can  be  computed  similarly  to  Eq.  11: 


EK\  =  al{l  +  {M-  1)[(1  +  02)2  +  M{  1  -  0!  -  02)[1  +  62  +  (1  -  -  02){2M  ~  l)/6]]}  (12) 

For  the  complete  model,  a  linear  frequency  drift  term  can  be  added.  That  is,  we  assume  that  one 
has  an  initial  data  set  (chosen  to  be  100  time  differences  between  a  pair  of  clocks).  The  number  of 
data  points  are  the  same  for  both  the  adaptive  servo  and  the  batch  processing.  The  introduction  of 
a  non-zero  frequency  drift  significantly  affects  both  the  adaptive  servo  and  the  conventional  Box  and 
Jenkins  analysis.  Equation  13,  below,  provides  the  additional,  independent  error  to  the  forecast  errors 
of  Eq.  12  due  to  frequency  drift: 

1  =  ffa{l  +  (1  -  #i)2  +  (#i  +  02)2  +  0\  +  (M  -  2)(1  —  61  -  02)2}/M2  (13) 

This  noise  addition  assumes  that  while  9\  and  02  are  known  exactly,  the  drift  term  is  estimated 
from  the  mean  second  difference  time  error  and  its  variance.  Appendix  A  contains  a  derivation  of  Eqs. 
12  and  13. 
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It  is  important  to  realize  that  the  actual  errors  in  real  applications  are  calculated  on  the  basis  of 
the  ESTIMATED  model  parameters,  not  the  “true”  parameters  which  are  unknown.  This  usage  of 
the  estimated  parameters  renders  the  computed  errors  a  bit  on  the  optimistic  side  because  Eq.  12 
does  not  include  the  errors  in  the  estimated  parameters.  Figure  7a  shows  the  theoretical  contributions 
of  the  conventional  Box  and  Jenkins  analysis  and  the  imperfect  knowledge  of  the  drift  rate  for  various 
theta  values.  Figure  7  gives  graphical  views  of  Eqs.  12  and  13. 

Figure  8a  shows  the  theoretical  forecast  errors  as  calculated  using  Eqs.  12  and  13.  For  simulation 
purposes,  the  “true”  parameter  values  are  known.  Simulation  techniques  allow  one  to  verify  Eqs.  12 
and  13  by  repeated  calculation  of  the  forecasts  compared  to  the  “true”  values  after  the  fact. 

The  bases  for  the  data  plotted  in  Fig.  8a,  b,  and  c  are: 

1.  ARIMA  (0,2,2)  model 

2.  Linear  frequency  drift 

3.  Data  length  of  100  points  (Initial  randomization  of  filter  required) 

4.  Forecast  100  lags  beyond  the  data  length  using  estimated  parameters 

5.  Repeat  above  for  independent  noise  samples  for  at  least  500  individual  runs. 

For  Fig.  8a,  the  errors  in  the  forecast  for  the  simulated  data  are  -.22  dB  worse  than  the  theoretical 
value.  (This  number  is  probably  within  the  uncertainty  limits  of  the  experiment.)  The  solid  line  is  the 
theoretical  result  of  Eqs.  12  and  13.  The  small  dots  near  the  solid  line  are  the  results  of  simulation. 

The  next  step  was  to  estimate  the  theta  parameters  using  a  conventional  Box  and  Jenkins  approach. 
Figure  8b  indicates  the  impact  on  the  forecasts  using  imperfect  parameters.  The  error  is  now  about 
1.01  dB  worse  than  optimum  after  100  lags.  Theoretical  errors  (i.e.,  solid  line)  are  identical  in  Figs. 
8a,  8b,  and  8c,  as  calculated  using  Eqs.  12  and  13. 

The  third  step  is  to  simulate  the  results  of  a  “real-time”  estimation  procedure  where  the  parameters 
are  “learned”  during  a  single  pass  through  the  100  data  points.  That  is,  the  adaptive  approach  is  used 
for  the  forecasts.  As  shown  on  Fig.  8c,  the  error  after  100  lags  is  now  1.14  dB  relative  to  optimal  or 
.13  dB  relative  to  the  Box  and  Jenkins  approach.  It  is  interesting  to  note  that,  at  least  for  the  given 
parameters,  the  adaptive  forecaster  is  close  to  the  Box  and  Jenkins  forecaster  in  performance. 

With  the  above  approach,  we  can  now  evaluate  the  adaptive  forecasting  performance  for  various 
physical  models  by  using  different  “true”  values  for  the  theta  parameters,  drift  rate,  and  noise  level. 
Further,  we  can  estimate  the  relative  performances  of  “perfect”  parameters,  Box  and  Jenkins  estima¬ 
tions,  and  “real-time”,  adaptive  method.  In  effect,  we  can  evaluate  the  costs  in  accuracy  of  using  the 
much  simpler  adaptive  approach. 


Reference 

•  G.E.R  Box  and  G.M.  Jenkins,  “Time  Series  Analysis”,  Holden-Day  San  Francisco,  1970. 
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APPENDIX  A 


PSI-WEIGHTS  FOR  AN  ARIMA  (0,  2,  1) 
Given  an  ARIMA(0,  2,  1)  Model: 


Xn  —  2  •  X„_i  +  Xn-2  —  an  +  0  ■  an-i  —  0 


A-l 


We  define  the  ^--weights  by  the  relation: 

OO 

Xn  =  ^2  \E\-  •  an- 1  for  n  =  0  to  oo 
i=o 


A-2 

In  order  to  have  both  (A-l)  and  (A-2)  valid  for  all  n,  we  substitute  (A-2)  into  (A-l)  and 
group  coefficients  of  the  a„  together  and  require  the  resulting  coefficient  to  vanish,  That  is, 
the  net  coefficient  of  an  is  just  'to  1.  The  first  few  relations  for  the  ^-weights  are: 

^o  =  l 
*1  =  (2-0) 

W2  =  (3-2-0) 

=  (4-3-0) 

'I'm  =  [rn  +  1)  —  m  ■  0 


A-3 


The  mean  square  forecast  errors  (See  Box  and  Jenkins)  are  given  by: 


M- 1 

■  m2  +  2  ■  (1  —  0)  •  m  +  lj  A-4 

TO= 0 

Ma\  U  +  (l-  0)(M  ~  1)  [1  +  (1  ~  0)(2M  -  l)/6]}  A-5 
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PSI-WEIGHTS  FOR  AN  ARIMA  (0,  2,  2) 


Given  an  ARIMA  (0,  2,  2): 


Xn  —  2X„_i  +  Xn~2  ~  0>n  +  Q\an~l  +  ^2°r»-2  —  0 


A-6 


We  define  the  ^-weights  as  is  ( A-2) .  The  first  few  relations  for  the  ^-weights  are: 


*1 

#2 

*3 


1 

2  -61 

2(2  —  #i)  —  (1  +  $2) 

3(2  -  $i)  -  2(1  +  02) 


__  (  m(  1  —  0i 

“  l1 


—  0i  —  0i)  +  (1  +  02)  for  m  >  0 

for  m  =  0 


The  mean  square  forecast  errors  become: 

E{6l,}  =  a\  {l  +  (M  -  1)  [(1  +  02)2  +  M(  1  -  0X  ~  02)  [1  +  02  +  (1  -  0i  -  02)(2M  -  l)/6]] } 
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APPENDIX  B 

Mean  and  variance  of  an  ARIMA  (0,  2,  2) 


Zn  =  a„  -  -  ^2«n-2  +  Aw 


B-l 


The  sum  of  the  first  M  values  of  Zn  are: 

M 

EZn  =  aM  +  (l  —  $l)  aM— 1  +  (1  ~  —  ^2)^-2  +  '  *  *  —  (#1  +  02)aO  “  #2a-l  +  M  •  Av 

n=l 


B-2 

Since  the  an  are  independent  random  numbers  with  zero  mean,  the  mean  value  is  obtained 
by  dividing  the  sum,  above,  by  M.  The  variance  of  the  estimate  is  obtained  by  taking  the 
expectation  value  of  the  square  of  (B-2).  That  is: 

V2  =  a2  {l  +  (1  -  0J2  +  (0!  +  02)2  +  (M  -  2)(1  -  61  -  02)2  +  02}  /M2 

B-3 

where  V2  is  the  variance  of  the  estimated  mean  of  Zn. 
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ARIMA  <0,  2,  2) 


FREQUENCY 


Fig.  1  Power  Spectra  f.or 


Parameter 


£s lima  ti 


on 


2b 


FIG. 4  ARIHA  FORECAST  SCHEME 
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ACTUAL  ARINA  COEFFICIENT,  THETA 


Fig.  5  Predominate  Errors  for  ARIMA(0,2,2) 
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7  FORECAST  TIME  ERRORS  FOR  AN  ARIMA  <0.  2.  2> 
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FORECAST  ERRORS  FOR  AN  ARIMA ( 0,  2 , 2  > 


RMS  ERROR 


RMS  ERROR 


RMS  ERROR 


Fig.  3a  Errors  Using  Exact  Parameters 


Fig.  8B  Errors  Using  Box/Jenkins  Estimation 


Fig  8C  Errors  Using  Adaptive  Estimation 
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QUESTIONS  AND  ANSWERS 

DR.  GERNOT  WINKLER,  USNO:  Since  you  can  do  your  parameter  estimation  on  the  run, 
you  can,  of  course,  allow  a  change  in  the  characteristics  of  your  frequency  standards  with 
time.  Its  an  adaptive  method. 

DR.  BARNES:  That  is  true.  It  will  adapt  to  the  new  value  if  the  standard  changes.  I  guess 
that  I  should  point  out  on  the  last  plot,  I  am  assuming  that  I  have  100  points  of  data  and 
I  am  forecasting  for  the  next  100  points.  It  is  not  a  lot  of  data  that  I  am  working  with. 
We  are  still  seeing  agreements  to  a  fraction  of  a  dB  in  this  approach. 
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